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INTRODUCTION 


The  special  vulnerability  of  the  head  to  serious  or  fatal  injury  is  well  documented  in  the  publication:  A 
Survey  of  Current  Head  Injury  Research”  prepared  by  the  Subcommittee  on  Head  Injury  of  the  Natior^ 
Advisory  Neurological  Diseases  and  Stroke  Council”"-  Broadly  speaking,  depending  on  the  state  of  the 
skull,  i.e..  whether  or  not  the  skull  is  penetrated  or  fractured  in  the  impact  process,  craniocerebral  trauma  is 
termed  either  open  or  closed.  This  investigation  is  concerned  with  the  injury  sustained  in  closed  head  impact. 

Models  of  the  closed  head  have  been  divided  in  much  the  same  fashion  as  the  major  hypotheses  proposed 
to  explain  closed  head  injuries  due  to  impact.  One  group,  the  rotational  school,  has  contended  that  the  rota- 
tional acceleration  induced  by  impact  causes  high  shear  strains  in  the  brain  matter,  rupturing  cerebral 
blood  vessels  and  tissue,  while  another  group,  the  cavitation  school,  has  claimed  that  there  exist  points 
within  the  brain  where  the  reduced  pressure  is  sufficient  to  rupture  the  capillary  walls.  Normally,  the  pres- 
sure differential  across  the  capillaries  is  only  a few  mmHg  During  the  impact,  the  pressure  fiela  at  a given 
time  contains  regions  of  negative  pressure,  i.e.,  less  than  atmospheric.  The  transmural  pressure  at  these 
locations  is  then  suddenly  increased,  which  could  lead  to  the  bursting  of  the  capillaries  therein.  Certainly, 
when  the  pressure  is  reduced  to  the  vapor  pressure  of  the  brain  substance,  cavitation  takes  place.  The 
catastrophic  collapse  of  the  bubbles  thus  formed  is  a possible  cause  of  brain  damage.  An  advanced  model 
of  the  cavitation  theory  was  given  by  this  author*”'-  Based  on  a sugge.stion  by  Goldsmith”"-  and  the 
ground  work  done  by  Anzclius”'-'- Guttinger*''''- ®',  Engin  and  Liu' modeled  the  head  as  a fluid-filled 
closed  spherical  shell  under  torsionless  axi.symmetric  impact  loading.  The  fluid  is  considered  inviscid  and 
represents  the  cerebrospinal  fluid  and  the  brain  matter.  The  shell  is  elastic,  homogeneous,  and  isotropic  and 
simulates  the  skull.  The  -ft'all  thickness  is  thin  but  includes  both  the  extensional  and  bending  effects.  The 
infinite-series  solution  was  obtained  by  Engin”-'  through  the  use  of  the  Laplace  transform  technique, 
for  the  case  of  an  impulsive  axisymmetric  polar-cap  load.  In  his  dissertation'”'*-  the  author  has  extended 
the  Engin  result  to  include  the  following  additional  considerations:  (a)  An  asymmetric  impact  loading 
consisting  of  an  axisymmetric  pulse  coincident  v/ith  a tangential  surface  torque  (b)  A moderately  thick  shell 
theory,  which  includes  extensional,  bending,  rotatory  inertia  and  transverse  shear  effects.  Engin  and 
Rob  arts” ” recently  examined  the  problem  of  the  transient  response  of  a fluid-filled  elastic  spherical 
shell  to  an  arbitrary  velocity  input  to  the  shell.  The  excess  pressure  distribution  in  the  fluid  was  evaluated 
for  various  deceleration  time  pulses.  As  a model  for  head  injury,  this  problem  belongs  to  the  noncontact 
variety.  Since  linear  deceleration  or  acceleration  of  the  shell  is  the  only  input,  the  mechanism  of  injury  can 
only  be  the  excess  pressure. 

As  for  models  of  rotation  theory,  Lee  and  Advani,  modeled  the  head  as  an  elastic  sphere  bonded  to  a rigid  shell 
undergoing  rotational  motion.  The  experimental  foundation  for  this  modelling  advance  is  typified  by  the 
recent  work  of  Ommaya  et  al.'^'-^’and  the  earlier  Holboum'”'*- Owings”'" working  on  work  of 
a dissertation  supervised  by  Prof.  S.  H.  Advani,  is  attacking  the  problem  of  the  dynamic  axisymmetric 
response  of  an  elastic  spherical  shell  with  an  elastic  core.  Hickling  and  \mer”"  modelled  the  head  as 
a two-layered  vi.scoelastic  sphere  using  the  exact  equations  of  linear  viscoelasticity.  The  loadings  and  motions 
are  considered  axisjonmetric.  The  above  problems,  because  of  their  formulation,  are  capable  of  evaluating 
the  relative  magnitude  and  duration  of  the  normal  and  shear  stresses  in  the  elastic  core.  When  the  material 
properties  and  failure  criterion  for  the  brain  and  skull  have  been  substituted  into  the  model,  a tolerance 
evaluation  can  be  made  for  a given  mechanism  of  injury. 

For  any  given  impact,  the  head  will  experience  both  a translational  and  rotational  acceleration.  Even  if  one 
were  to  control  the  blow  in  such  a manner  that  the  force  vector  passes  through  either  the  center  of  mass 
or  center  of  percussion  of  the  head,  both  motions  will  occur  because  of  the  constraint  exerted  by  the  head- 
neck  junction.  Implicit  in  many  of  the  papers  of  both  the  translational  and  rotational  theory  of  brain  damage 
is  the  impression  that  translational  motions  induce  only  nonnal  stresses  or  pressure  and  in  rotational  motion, 
only  shear  stresses.  This  is  not  generally  true.  From  continuum  mechanics  it  is  known  that  a normally 
applied  load  induces  shear  stresses  in  all  but  the  normal  directions.  Conversely,  tangential  surface  tractions 
induce  normal  stresses  using  a similar  argument.  'Which  effect  is  primary  during  a given  impact  is  the  central 
question. 
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It  is  tho  intent  of  this  invcstifjntion  to  dotornninc  the  relative  magnitude  of  these  two  effects  under  a given 
impact.  To  do  so,  a more  realistic  model  has  to  he  ctmeeived.  Tho  inviscid  fluid  a.ssumt'tion  mu.st  give  way 
to  include  viscosity  and  a multi-layered  ellipsoiflal  shell  should  replace  the  spherical  shell.  In  so  doing,  the 
complexities  of  the  boundary  conditions  together  with  the  nonlinearities  in  the  enclosed  flow  make  a 
closed-form  solution  nearly  impossible.  The  approach,  which  shows  the  most  promise  for  this  analytically 
intractable  problem,  is  the  finite  element  method. 
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FORMULATION  OF  THE  PROBLEM 


Finite  Element  Method  The  triuiitional  method  of  attacking  a physical  problem  consists  jn  choosing  a 
suitable  “infinitesimnl”  element  in  the  continuum  and  foimulating  a differential  or  integral  equation  for 
it.  Such  an  equation,  though  derived  from  physical  laws,  is  not  itself  a physical  law;  it  is  a formulation, 
derived  by  a certain  amount  of  prior  mathematical  deduction.  The  finite  element  method  is  a relatively 
new  technique  in  physical  science.  It  is  flistinctly  different  from  conventional  mathematical  or  numerical 
analysis.  We  cannot  say  that  this  type  of  analysis  is  inherently  better  than  any  other  as  it  does  have  some 
limitations,  but  it  also  has  some  distinct  advantages. 

Finite  element  analysis  suggests  that  we  move  the  power  of  numerical  analysis  forward  in  the  scheme  of 
things,  giving  it  a more  direct  or  even  central  role  in  the  analytic  process.  It  is  essentially  a process  through 
which  a continuum  system  with  infinite  degrees  of  freedom  can  be  approximated  to  by  an  assemblage  of 
.subregions  or  elements  each  with  a specified  but  finite  number  of  unknowns.  In  this  study  the  displacement 
formulation  of  the  finite  element  method  has  been  applied.  We  start  by  dividing  the  continuum  geometrically, 
in  some  convenient  manner,  into  a finite  number  of  elements  of  appropriate  size.  The  elements  are  assumed 
to  be  interconnected  at  a di.screte  number  of  nodal  points  situated  on  their  boundaries.  The  displacements 
of  these  nodal  points  are  the  basic  unknown  variables  of  the  problem.  A set  of  functions  called  displacement 
functions,  is  chosen  to  define  uniquely  the  state  of  displacement  within  each  finite  element  in  terms  of  its 
nodal  displacem.ent.  Each  of  these  elements  is  then  expre.ssed  in  such  a manner  that  it  satisfies  the  basic 
laws  and  the  kinematic  and  constitutive  relations,  together  with  certain  assumptions  and  boundary  condi- 
tions. WTiatever  changes  a system  may  undergo  in  a certain  interval  of  time,  this  method  of  analysis  demands 
the  straight  forwa  rd  criterion  that  at  each  instant  every  individual  element  must  satisfy  a balance  of  mass, 
momentum,  and  energy,  and  must  bear  pertinent  kinematic  relationship  to  its  neighboring  elements.  In 
fact  the  finite  element  method  is  equivalent  to  the  minimization  of  the  total  potential  energy  of  the  system 
in  terms  of  a prescribed  displacement  field.  Its  application  can  be  extended  to  those  problems  where  a 
variational  formulation  is  possible. 

Material  Properties  and  Constitutive  Equations  The  human  head  is  modeled  as  a solid  viscoelastic  core 
bonded  lO  a viscoelastic  spherical  shell,  which  .simulates  the  brain  mattei  and  the  skull,  respectively.  Linear 
viscoelastic  properties  are  assumed  for  both  the  brain  matter  and  the  skull  along  with  homogeneity  and 
isotropy  assumptions.  In  terms  of  stress  matrix jr,  strain  matrix j_.  and  strain  rate  matrix the  constitutive 
equations  for  the  brain  matter  and  the  skull  for  the  axisymmetric  case  in  polar  coordinate  system  are  in  the 
form: 


= D I r + D:.  t 
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D,  and  D , are  4x4  matrice.s,  which  are  functions  of  the  material  properly  constants. 
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Derivation  of  The  Equations  of  Motion  In  the  displacement  formulation  of  the  finite  clement  method, 
the  basic  unknown  variable  are  the  displacements.  The  state  of  displacement  u in  a continuum  system 
is  approximated  in  terms  of  the  nodal  displacements  U by  a set  of  chosen  displacement  functions_f,  i.e. 


u = fU 


(2) 


3 


where 


rut 

r 


u = 


, U = 


u, 

u., 


u,. 


and 


h,(r,z)  h.:(r.z)  • • • h,.fr,z) 

f(r,z)=  |k,(r,z)  kj(r,z)  • • • k„(r,z) 
Applying  the  principle  of  virtual  wcrk  for  dynamic  systems,  i.e. 

SE  = SW  — u dv 


(3) 


where  E and  W are  the  virtual  strain  energy  and  virtual  work,  respectively,  for  any  compatible  virtual 
displacement  The  last  term  in  the  above  equation  is  the  work  due  to  inertial  force,  and  /^is  the  density 
of  the  material. 

Taking  first  and  second  derivatives  with  respect  to  the  independent  variable  time  and  spatial  variables  of 
equation  2,  we  have 

u = fij  , u = fU 

and  (4) 

= gU  , £ = glj_ 

where  g = g(r,z)  is  obtained  by  taking  first  derivatives  with  respect  to  the  independent  spatial  variables 
of  the  displacement  function  matrix  f(r,z),  specifically. 
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From  equations  1 and  4,  we  obtain 
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(t  = D,  gU  + D gU 

Neglecting  the  body  force,  then  the  virtual  strain  energy  and  virtual  work  due  to  surface  load  p,  are 
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B = I tVdv 

and  (6) 

W = / u^'p  ds 


Substituting  equation  6 into  equation  3 and  making  use  of  equations  4 and  5,  we  have 

( [5U]^  g'-[D,gU  + D.gU]  dv 
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Since  SU  is  arbitrary. 


[/  gj  D,  g dv]  U + [/  _gj  D=ldv]  U 
=7ppds-[/Fp£dv]  U 

s V 

MU  + CU  + KU  = P (7) 

where 

M — J F p f dv  represents  the  mass  matrix  of  the  equivalent  discrew;  system 
C=j  g'*' Dl. g dv  is  the  viscous  damping  matrix 
K=/  g''  D,gdv  is  the  stiffness  matrix,  and 

V ~ ~ 

F =J  P pds  is  the  equivalent  concentrated  force  due  to  the  surface  load. 


METHOD  OF  S0LUTI0^4 

The  governing  equations  of  the  dynamic  system  are  expressed  by  the  second  order  differential  equations 
in  matrix  form,  shown  in  equation  7,  which  are  essentially  a state  of  equilibrium  between  the  internal  and 
the  external  forces  on  each  element  of  the  system.  Solutions  to  these  equations  are  accomplished  through 
an  iterative  procedure  to  generate  the  time  histories  by  first  determining  the  insta  itaneous  acceleration 
over  a short  time  interval  for  each  degree  of  freedom  and  then  perfoiming  integrations  to  obtain  the  velocity 
and  displacement  changes  over  the  short  time  interval.  These  velocity  and  displacement  changes  are  then 
added  to  their  previous  values  to  update  the  time  histories.  The  brief  outline  for  the  procedure  is  as  follows: 
Assuming  the  acceleration  is  linear  over  the  small  time  increment  At  = t„  — tn.,,  then 


(i)  The  velocity  at  t„  is  calculated  as 

The  computation  starts  with  an  assumed  value  of  = 0 (or  any  other  reasonable  values) 

(ii)  The  displacement  at  t„  is  then  computed  from 

= y.„„  + 1/2  (At)  + y. j + i/fe  (-..)--{b7.  - Uu. ,} 

• • 

(iii)  The  error  introduced  by  the  assumed  values  of  U,,,  in  step  (i)  is  coro;cted  by  re-onmputirq  ;;  v accelera- 
tion y,„  from  the  equations  of  motion,  eq.  7,  with  the  velocity  U-.  and  rlisplactmeut  U,.  giver;  bv  (i)  and 
(ii)  i.e. 

My,.zzP.„-cy.,.-Ky,„ 

• • 

with  this  computed  value _U,,^  the  calculations  go  back  to  steps  (i)  through  (iii)  and  thi?  pr(>ceduro  is 
repeated  again  until  convergence  has  been  obtained  to  the  desired  accuracy. 


NUMERICAL  RESULTS  AND  OBSERVATIONS 

The  geometric  dimensions  and  the  material  constants  used  for  the  inner  solid  viscoelastic  core,  which  repre- 
sents the  brain  matter,  and  the  outer  viscoelastic  spherical  shell,  which  simulates  the  skull,  are  as 
follows'""'- 
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density 
bulk  modulus 
shear  modulus 
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the  applied  rec^an^nilar  pulse  Is  easily  observed.  Both  tb.e  ampiitudf  and  the  pattern  of  the  pressure  (stress) 
field  in  the  core  are  altered  wuh  the  ir.creascd  diiraticri.  The  pressr.re  histories  at  the  impact  pole  for  pulses 
of  different  dui'-tiori.-  are  (-nvo’',  ir-  ftyurrs  7 ;-md  f’.  From  these  and  tne  foregoing  figures,  it  is  evident  that  the 
increase  in  pulse  length  not  oaly  moreases  the  amplitude  of  the  pre.vsure  field  but  also  increases  the  duration 
of  the  rarefaction,  i.e..  reduce  ' pressure  or  tensile  stress,  in  the  core.  Rarefaction  or  reduced  pressure  is  the 
cause  of  brain  damage  according  to  the  cavitation  hypothesis.  The  apparent  wave  speed  in  the  viscoelastic 
core  can  be  calculated  from  figures  5 or  6 by  observing  the  arrival  time  of  the  wave  front  at  different  locations 
and  the  distance  between  these  locations.  The  calculated  wave  speed  is  53,700  in./sec  which  is  very  close  to 
the  esperimental  wave  speed  in  brain  matter"’"  ”*. 

In  general,  the  duration  of  an  actual  im.pact  on  the  head  is  approximately  one  millisecond  or  longer  and  the 
two  rectangular  pulses  with  durations  of  .01  and  .06  msec  respectively,  as  used  in  figures  7 and  8 are  not 
realistic.  They  are  used  here  only  to  demonstrate  the  effect  of  the  impact  duration  on  the  .stress  (or  pressure) 
field,  rather  than  to  simulate  a physically  obtainable  impact.  A more  meaningful  and  realistic  situation  is 
simulated  by  a pulse  with  an  exponential  time  function  as  shown  in  figure  4.  Results  obtained  with  this 
loading  pulse  are  presented  in  figures  9 through  15. 

The  pressure  histories  of  locations  where  rarefaction,  i.e.,  reduced  or  negative  pressure,  occur  are  given  in 
figures  9 and  10.  All  those  locations  are  situated  along  the  axis  of  symmetry,  or  the  impact  diameter  (The 
diameter  passing  through  the  centroid  of  the  impact  area).  These  correspond  to  possible  locations  of  brain 
damage  according  to  the  cavitation  hypothesis.  The  pressure  histories  for  the  locations  at  2=3.0  in.  (the 
impact  pole),  z=0.12  in,  (close  to  the  center  of  the  sphere),  z=  _ 1.67  in.,  and  z=  —3.0 in.  (the opposite 
pole),  are  given  in  figures  9 through  12,  respectively.  From  these  figures,  it  can  be  seen  that  the  maximum 
rarefaction  occurs  at  the  counter  pole,  which  is  followed  by  the  one  situated  at  2=  — 1.67  in.  This  one  is 
moderate  in  magnitude  and  last.'?  longer.  Figure  12  shows  that  the  opposite  pole  is  under  negative  pressure 
(tension)  during  a great  part  of  the  impact  period. 

As  for  shear  .stre.sses  in  tlie  viscoelastic  core,  from  the  exi>onential  time  pulse,  there  are  two  regions  where 
the  shear  stresses  are  significant.  One  is  located  along  a line,  approximately  parallel  to  and  0.4  in.  away 
from  the  axis  of  symmetry.  The  otlier  one  is  situated  along  the  circumference  of  a 2.2  in.  radius  circle,  i.e.. 
approximately  0.8  in,  from  e .surface  of  the  shell.  The  shear  stresses  in  these  two  regions  are  given  in 
figures  13  and  14  respectively.  In  figure  13  the  four  curves  represent  the  shear  stresses  at  times  t = 0.16, 
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0.56,  2.48,  and  3.34  n>sec,  rcspnctivoly,  While  ihc  three  eun'c-;  h,  fipure  M !;ho\^  I he  Khear  ptrfc.->s  at  tiroes 
t=:0.18,  0.62,  and  1.42  msec,  respectively.  The  maximum  sbtsar  stresM  i.s  .<iii]nv;HijmateIy  1.2  lb/in.‘.  With 
shear  stress  of  this  magnitude,  these  t\vt>  regions  are  pos.sible  locations  oJ  Inain  damage,  according  to  the 
rational  hypothesis. 

The  stress  distributions  in  the  shell  are  shown  in  Fig'jre  15.  The  fhrei;  cur.^es  thet?  represent  the  normal 
stress  in  the  polar  angle  direction  at  times  t n;  O.iu,  i.68,  and  2.72  inst'c,  respt:et:ve}y.  From  these  ii  ivui  be 
seen  that  stresses  at  the  tw'o  poles  are  relatively  large.  It  i.-j  gcner/illy  l>eHev‘ d <hai  skull  fractures  are 
due  to  high  tensile  stresses;  and  in  vnew  of  this,  a .skull  fracture  would  ()r.7b.alijy  rtreur  ai  either  pole  or  in 
the  neighborhood  of  a cone  with  an  exiended  angle  of  50'.  Further,  both  l>ca.iion.s  also  have  high  sherir 
stresses. 
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Figure  1.  Finite  Element  Mesh 
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of  Symmetry  for  Square  Pulse  of  0.06  msec  Duration 
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Figure  7.  Pressure  of  Brain  at  the  Impact  Pole  for  Square  Pulse  of  .01  msec  Duration 


Figure  8.  Pressure  of  Brain  at  the  Impact  Pole  for  Square  Pulse  of  0.06  msec  Duration 
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Fif'ure  l.'l.  Shear  Stress  History  Along  the  Diametric  Axis  for  the  Exponential  Time  Input  Pulse 
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